In [10]:
from __future__ import division, print_function
# Standard library
import os
# Third-party
from astropy import log as logger
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
# Custom
import streamteam.dynamics as sd
import streamteam.integrate as si
import streamteam.io as io
import streamteam.potential as sp
from streamteam.units import galactic
np.set_printoptions(linewidth=128)
In [34]:
def read_aaf(path, snapfile):
actions = np.load(os.path.join(path,"actions_{}.npy".format(snapfile)))
angles = np.load(os.path.join(path,"angles_{}.npy".format(snapfile)))
freqs = np.load(os.path.join(path,"freqs_{}.npy".format(snapfile)))
return actions,angles,freqs
In [35]:
actions,angles,freqs = read_aaf("/Users/adrian/projects/morphology/simulations/runs/thin3/actions/", "SNAP060")
In [36]:
dJ = actions[1:] - actions[0]
dF = freqs[1:] - freqs[0]
In [37]:
dF_mag = np.linalg.norm(dF, axis=1)
leading = np.all(np.sign(dF) < 0., axis=1)
In [13]:
def find_hessian_eigvals(dF, dJ):
N = dJ.shape[0]
Y = np.ravel(dF)
A = np.zeros((3*N,9))
A[::3,:3] = dJ
A[1::3,3:6] = dJ
A[2::3,6:9] = dJ
# Solve for 'parameters' - the Hessian elements
X,res,rank,s = np.linalg.lstsq(A, Y)
# Symmetrize
D0 = X.reshape(3,3)
D0[0,1] = D0[1,0] = (D0[0,1] + D0[1,0])/2.
D0[0,2] = D0[2,0] = (D0[0,2] + D0[2,0])/2.
D0[1,2] = D0[2,1] = (D0[1,2] + D0[2,1])/2.
print("Residual: " + str(res[0]))
return D0,np.linalg.eigh(D0) # symmetric matrix
In [38]:
l_hess,(l_eigvals,l_eigvecs) = find_hessian_eigvals(dF[leading], dJ[leading])
t_hess,(t_eigvals,t_eigvecs) = find_hessian_eigvals(dF[~leading], dJ[~leading])
In [39]:
print(l_eigvecs)
print(l_eigvals)
In [40]:
print(t_eigvecs)
print(t_eigvals)
In [41]:
e1,e2,e3 = l_eigvecs.T
Q,R = np.linalg.qr(np.outer(e1, e1)/l_eigvals[0] + np.outer(e2, e2)/l_eigvals[1] + np.outer(e3, e3)/l_eigvals[2])
print(np.linalg.inv(R).dot(Q.T))
print(l_hess)
In [42]:
print(np.linalg.inv(l_hess))
print(np.outer(e1, e1)/l_eigvals[0] + np.outer(e2, e2)/l_eigvals[1] + np.outer(e3, e3)/l_eigvals[2])
In [43]:
dims = (0,2)
for i in range(3):
ei = l_eigvecs[:,i]
plt.plot([0,ei[dims[0]]], [0,ei[dims[1]]], linewidth=2, marker=None)
Another idea: integrate a "ball" of orbits around the initial orbit of interest, use those to compute Hessian!
In [179]:
def orbit_ball(w0, N, scale=0.01):
""" scale sets the relative size of the ball (0.01 = 1%) """
R = np.linalg.norm(w0[:3])
V = np.linalg.norm(w0[3:])
scale = np.array([R,R,R,V,V,V]) * scale
w0s = np.random.normal(w0, scale, size=(N,6))
return np.vstack((w0,w0s))
# w0 = [20.,2.5,4,0.,0.1,0.17]
# grid = orbit_ball(w0, N=100, scale=0.007)
pal5_w0 = [8.312877511, 0.242593717, 16.811943627] + list(([-52.429087, -96.697363, -8.156130]*u.km/u.s).to(u.kpc/u.Myr).value)
grid = orbit_ball(pal5_w0, N=256, scale=1E-3)
In [180]:
ppars = {'a': 6.5,
'b': 0.26,
'c': 0.3,
'm_disk': 65000000000.0,
'm_spher': 20000000000.0,
'phi': 1.570796,
'psi': 1.570796,
'q1': 1.3,
'q2': 1.0,
'q3': 0.8,
'r_h': 30.0,
'theta': 1.570796,
'v_h': 0.5600371815834104}
# potential = sp.PW14Potential(**ppars)
potential = sp.LM10Potential()
In [181]:
t,w = potential.integrate_orbit(grid, dt=0.4, nsteps=70000, Integrator=si.DOPRI853Integrator)
In [182]:
fig = sd.plot_orbits(w, alpha=0.01, linestyle='none')
In [183]:
actions,angles,freqs = sd.find_actions(t[::20], w[::20], N_max=8, units=galactic)
In [199]:
fig,axes = plt.subplots(1,3,figsize=(15,5))
q = freqs*1000.
q -= q.min(axis=0)
axes[0].plot(q[:,0], q[:,1], linestyle='none')
axes[1].plot(q[:,0], q[:,2], linestyle='none')
axes[2].plot(q[:,1], q[:,2], linestyle='none')
fig.suptitle(r"Projections of $\Omega - \mathrm{min}(\Omega)$ [Gyr$^{-1}$]", fontsize=24, y=1.02)
Out[199]:
In [185]:
fig,axes = plt.subplots(1,3,figsize=(15,5),sharex=True,sharey=True)
q = (angles % (2*np.pi)) / (2*np.pi)
axes[0].plot(q[:,0], q[:,1], linestyle='none')
axes[1].plot(q[:,0], q[:,2], linestyle='none')
axes[2].plot(q[:,1], q[:,2], linestyle='none')
axes[0].set_xlim(0,1)
axes[0].set_ylim(0,1)
Out[185]:
In [150]:
hess,(eigvals,eigvecs) = find_hessian_eigvals(freqs[1:]-freqs[0], actions[1:] - actions[0])
In [151]:
hess
Out[151]:
In [152]:
eigvals
Out[152]:
Mean fundamental frequency is of order 1.4E-2 Myr^-1, so mean orbital period is about 70 Myr.
I'll integrate an orbit for 2 x 100 orbital times, split it in two and FFT both
In [163]:
Tint = 70 * 2000 # Myr
dt = 1.
# fine_t,fine_w = potential.integrate_orbit(w[0,0], nsteps=int(Tint/dt), dt=dt, Integrator=si.DOPRI853Integrator)
fine_t,fine_w = t,w
In [165]:
fig = sd.plot_orbits(fine_w, ix=0, linestyle='none', alpha=0.1)
In [166]:
nsteps = fine_w.shape[0]
split = nsteps//2
f_1,fft_1 = sd.nonlinear.fft_orbit(fine_t[:split], fine_w[:split,0])
f_2,fft_2 = sd.nonlinear.fft_orbit(fine_t[split:], fine_w[split:,0])
In [167]:
f_2[:len(f_2)//2]
Out[167]:
In [170]:
plt.figure(figsize=(12,12))
plt.loglog(f_1[:,0], fft_1[:,0], alpha=0.7)
plt.loglog(f_2[:,0], fft_2[:,0] / fft_2[len(f_2)//2,0] * fft_1[len(f_1)//2,0], alpha=0.7)
# plt.xlim(1E-4, 1E-2)
# plt.ylim(4E2, 10**(5.5))
plt.title("FFT")
for freq in freqs[0]:
plt.axvline(freq, color='r') # fundamental frequencies
In [169]:
plt.figure(figsize=(12,12))
plt.semilogx(f_1[:,0], fft_1[:,0]**2, alpha=0.7)
plt.semilogx(f_2[:,0], fft_2[:,0]**2, alpha=0.7)
# plt.xlim(-1E-1, 1E-1)
# plt.ylim(4E2, 10**(5.5))
plt.title("Power spectrum")
for freq in freqs[0]:
print(freq)
plt.axvline(np.abs(freq), color='r') # fundamental frequencies
In [ ]: